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Abstract. The isotropic-to-nematic transition in a two-dimensional fluid of hard needles is studied using 
grand canonical Monte Carlo simulations, multiple histogram reweighting, and finite size scaling. The 
transition is shown to be of the Kosterlitz-Thouless type, via a direct measurement of the critical exponents 
Tj and /3, of the susceptibility and order parameter, respectively. At the transition, = 1/4 and /3 = 1/8 are 
observed, in excellent agreement with Kosterlitz-Thouless theory. Also the shift in the chemical potential of 
the nematic susceptibility maximum with system size is in good agreement with theoretical expectations. 
Some evidence of singular behavior in the density fiuctuations is observed, but no divergence, consistent 
with a negative specific heat critical exponent. At the transition, a scaling analysis assuming a conventional 
critical point also gives reasonable results. However, the apparent critical exponent P^s obtained in this 
case is not consistent with theoretical predictions. 

PACS. 64.60.Pr Equilibrium properties near critical points, critical exponents - 64.70.Md Transitions 
in liquid crystals - 64.60.Cn Order-disorder transformations; statistical mechanics of model systems - 
05.50.+q Lattice theory and statistics (Ising, Potts, etc.) 



1 Introduction 

Upon increasing density, a fluid of hard needles in two 
dimensions undergoes a transition from an isotropic to a 
nematic phase [T|. In the isotropic phase, the orientational 
correlations decay exponentially to zero, while in the ne- 
matic phase algebraic decay is observed Q- In the thermo- 
dynamic limit, long-range nematic order is thus absent in 
both phases, and the available evidence points to a tran- 
sition of the Kosterlitz-Thouless (KT) type |1|2] . In other 
words, the universality class of the isotropic-to-nematic 
(IN) transition in two-dimensional hard needles should be 
that of the XY model [3], and one expects to find the 
same set of critical exponents. For the XY model, the lat- 
ter are known exactly, but their verification in a fiuid of 
hard needles remains elusive to this day. The purpose of 
this paper is to fill this gap, using grand canonical Monte 
Carlo simulations and finite-size scaling. Indeed, our sim- 
ulations consistently recover the XY exponents 77 = 1/4 
and (3 = 1/8, of the susceptibility and order parameter, 
respectively. In addition, the scaling of the chemical poten- 
tial at the susceptibility maximum is in good agreement 



^ For this reason, the nematic phase should perhaps be 
termed quasi-nematic. For notational convenience, however, we 
refrain from doing so in this paper. 



with XY universality. Hence, our data quantitatively con- 
firm the KT scenario in fluids of hard needles. 

We also observe that, at high density in the nematic 
phase, the decay of nematic order with increasing sys- 
tem size is very slow. This means that even in macro- 
scopic samples a substantial degree of nematic order is 
present. The same occurs in the XY model: even though 
long-range magnetic order is absent in the thermodynamic 
limit, finite XY systems at low temperature nevertheless 
reveal considerable magnetic order. The consequences of 
this have been worked out by Bramwell and Holdsworth 
(BH) , who conclude that the formation of magnetic order 
in finite XY models is characterized by an effective critical 
exponent (3cS [Ij. Interestingly, when we analyze our data 
assuming a conventional critical point, we can also consis- 
tently measure such effective exponents, which moreover 
obey the hyperscahng relation. Hence, the BH scenario for 
the XY model seems to be valid in fluids of hard needles 
also, even though /3cff obtained by us differs from the XY 
value predicted by BH. 

The outline of this paper is as follows. We flrst specify 
the model details and the simulation method. Next, we 
present our raw simulation data, displaying how the var- 
ious observables of interest depend on the chemical po- 
tential and system size. The raw data is then analyzed 
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using several finite size scaling methods. We end with a 
discussion and summary. 



2 Model and simulation method 

We consider a two-dimensional system of infinitely thin 
rods of length I, henceforth referred to as needles. We em- 
phasize that our model is not discretized in any way: both 
the needle positions and orientations are continuous. In 
what follows, I will be the unit of length. The needles are 
hard, i.e. they are not allowed to overlap, and trivial fac- 
tors of inverse temperature are set to unity throughout. 
The simulations are performed in the grand canonical en- 
semble, i.e. at constant chemical potential ^ and system 
area A, while the number of needles N fiuctuates. The 
average needle density increases with ^ and this can be 
used to induce the IN transition. Hence, ^ is the con- 
trol parameter, analogous to inverse temperature in ther- 
motropic systems. We use a two-dimensional simulation 
square of size A = L"^ with periodic boundary conditions. 
Insertion and removal of needles are attempted with equal 
probability, and accepted with the standard grand canon- 
ical MetropoHs probabilities |Sj6j- During insertion, a ran- 
dom location in the system is selected and a needle with 
randomly selected orientation is tentatively placed at this 
location. If this needle overlaps with any of the other nee- 
dles already present, the move is rejected. Otherwise, the 
new state is accepted with probability 



A{N N +1) 



1, 



N + 1 



(1) 



with N being the number of needles in the system at the 
start of the move. Similarly, during removal, one of the 
needles is selected at random and deleted from the system, 
and the resulting state is accepted with probability 



A{N N -1) 



1, 



N 
Ae^' 



(2) 



To facilitate the efficient detection of overlap during parti- 
cle insertion, a link-cell neighbor Hst is used [7]. The sim- 
ulation data are collected as two-dimensional histograms 
Hfi^LiS, N), counting how often a state with nematic or- 
der parameter S and particle number N is observed (note 
the dependence on ^ and L). For system sizes L = 10 — 30, 
histograms are obtained for several values of fi; the mul- 
tiple histogram method [8|9j is used to evaluate proper- 
ties at intermediate values. The nematic order parameter 
S is defined as the maximum eigenvalue of the orienta- 
tional tensor Qaf3 — X^ili {'2'diadip — Sap)-, with dia the a 
component (a = x,y) of the orientation of molecule i, 
\di \ = 1, and 5ap the Kronecker delta. We emphasize that 
the nematic order parameter S defined in this way is an 
extensive quantity. In cases where the number of particles 
N is constant, it is convenient to use the normalized inten- 
sive definition S* = S/N ^ since then one has S* = Q and 
S** = 1, in an isotropic and perfectly aligned sample, re- 
spectively. However, in the grand canonical ensemble, N is 



a fiuctuating quantity, which itself might exhibit singular 
behavior, and so this convention is not used here. 

Most of the simulations were performed on Intel Dual- 
Core processors clocked at 2 GHz. For each system size L, 
around 10 histograms were collected, with ^ taken from 
the range 5.0 — 5.2. At these values, the needle density 
p ^ 7, which is close to the transition density observed 
in previous studies |1I2) . Each histogram was simulated 
for « 10^ grand canonical sweeps [10], with a sweep being 
defined as one complete renewal of the particle popula- 
tion (recall that the number of particles fluctuates). The 
computational effort per sweep depends on /i and L. For 
fx = 5.1 and L = 10, n « 2.3 x 10^ grand canonical Monte 
Carlo attempts are required to complete one sweep; for 
L = 30, this increases to n « 2.7 x 10^. The simulations 
began with empty boxes, and the first 1000 sweeps were 
discarded for equilibration. 



3 Results 
3.1 Observables 

The observables of interest are the average needle density 
and the compressibility 

p=(iV)M, Xp={{N^)-{Nf)lA, (3) 

the nematic density (order parameter) and the nematic 
susceptibility 

a = (5)M, X. = {{S^) - {Sf) I A, (4) 
and the Binder cumulant 



C/4 = (5^)7(^*> 



(5) 



The above quantities will generally depend on /i and L, 
especially in the vicinity of phase transitions. 



3.2 Raw simulation data 

We first present our raw simulation data. In Fig. [1] the 
nematic susceptibility x<t is plotted versus /i for several 
system sizes. The nematic susceptibility displays a maxi- 
mum, which becomes more pronounced as L increases. In 
addition, the chemical potential at the maximum /x^ de- 
pends on L. The nematic density is shown in Fig. [21 We 
observe that a increases monotonically with ^, while it 
decreases with increasing L. The Binder cumulant C/4 is 
shown in Fig. [3l The data from the different system sizes 
approximately intersect. The compressibility Xp is shown 
in Fig. [31 As with Xo-i the formation of a maximum is 
visible, but it increases only mildly with L. In Fig. [H the 
needle density is shown. We observe a monotonic increase 
of p with /i, and a weak decrease of p with L. In Fig. [H 
we show the nematic density at high values of the chem- 
ical potential, chosen well beyond the extrema of Xu and 
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5.20 



Fig. 1. Plot of the nematic susceptibility Xct versus the chem- 
ical potential /i, for several system sizes L as indicated. Clearly 
visible is that Xa attains a maximum, and that the maximum 
grows with increasing system size. Note also that the position 
of the maximum is size dependent. 




0.915 



5.20 



Fig. 3. Plot of the Binder cumulant Ui versus the chemical 
potential /i, for several system sizes L as indicated. Note that 
the data from the various system sizes approximately intersect. 
For increasing L, a shift of the intersection point toward higher 
values of Ui is visible. 




5.20 



Fig. 2. Plot of the nematic density a, which plays the role 
of the order parameter, versus the chemical potential /i, for 
several system sizes L as indicated. Note that a increases with 
/i, but also that it overall decreases with increasing L. This 
result is compatible with the absence of nematic order in the 
thermodynamic limit. 
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Fig. 4. Plot of the density fluctuation (compressibility) Xp 
versus the chemical potential /x, for several system sizes L as 
indicated. Note the presence of the maximum, but also that 
the increase of the maximum with L is much milder compared 
to that of X" in Fig- E 



Xp. We observe that a does not saturate, but continues to 
decrease with increasing L. 

The raw simulation data already provide evidence of a 
phase transition. Based on Fig. [H the transition is char- 
acterized by a diverging nematic susceptibility. There is 
also evidence of singular behavior in the compressibility, 
see Fig. IH but it is much weaker. The key observation is 
that the transition does not yield any finite order param- 
eter: even at very high chemical potential, the nematic 
density a does not saturate, but continues to decay with 
increasing L, see Fig. [6l This is consistent with previous 
simulations of hard needles |1|2) , and provides further con- 
firmation that nematic order is most likely absent in the 



thermodynamic limit, i.e. a decays to zero as L ^ c« 
irrespective of /z. The absence of nematic order appears 
to be a general property of two-dimensional liquid crys- 
tals - computer simulations of rods reveal similar behavior 
|ll|12j - and is conform the Mermin- Wagner theorem |13j . 
Note that for certain liquid crystal pair potentials, the ab- 
sence of nematic order can be proved rigorously [Hj. The 
observation that the order parameter vanishes in the ther- 
modynamic limit rules out a conventional critical point, 
leaving a transition of the KT type as the most likely al- 
ternative. 
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Fig. 5. Plot of the particle density p versus the chemical 
potential for several system sizes L as indicated. The overall 
trend is that p increases with /i, and that it decreases mildly 
with L. 




0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 
1/L 

Fig. 6. Variation of the nematic density a (order parameter) 
versus 1/L measured at two high values of the chemical poten- 
tial p, as indicated (by high is meant well beyond the maxima 
in Xct and Xp)- The important result to take from this graph 
is that a does not saturate at a finite value, but continues to 
decrease with increasing L, consistent with the absence of true 
nematic order in the thermodynamic limit. 
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Fig. 7. Determination of the thermodynamic limit transi- 
tion chemical potential pao assuming the KT scenario. Plot- 
ted is the chemical potential p*j^ of the nematic susceptibility 
maximum versus l/(lnL)'^; the line is a fit to the KT form of 
Ea. l|10p . The data follow the KT prediction well, and from the 
fit ^ 5.187 is obtained. 
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Fig. 8. Determination of the thermodynamic limit transition 
density poo assuming the KT scenario. Plotted is the density 
pL obtained At p = p^ versus l/(lnL)^; the line is a fit to 
Ea. lflO)l . Prom the fit pcx> ~ 6.98 is obtained. 



4 Finite size scaling 
4.1 KT scaling 

Characteristic of a KT transition is the exponential di- 
vergence of the correlation length fSj. If one starts in the 
isotropic phase, and moves toward the nematic phase by 
increasing the chemical potential, the correlation length 
diverges as 

^ oc exp ,t = Poo - fJ.,t>0, (6) 

with Poo the chemical potential at the transition, and 
nonuniversal constant 6 > 0. As ^ diverges faster than 



any power law, the conventional critical exponent v of the 
correlation length does not exist, but it is possible to de- 
fine exponents f3 and rj, of the nematic order parameter 
a and susceptibility Xtr, respectively, by expressing these 
quantities in terms of ^ 

cTcxr^ x<x«f'-^ (7) 

with KT values /? = 1/8 and rj = 1/4 [5]. We emphasize 
that these exponents are only observed on the positive 
interval < i < e, with e not too large [I 5;l6j . In the 
regime t < 0, the correlation length remains infinite, and 
the exponents become functions of t. 

In finite systems, the L dependence of a and Xa- in the 
regime < t < e is described in the context of finite size 
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Fig. 9. Application of Loison's method [20] to obtain the 
critical exponent in a two-dimensional fluid of hard needles. 
Plotted is Xo-(-f') Z/~'^~''' versus U4, for various system sizes L, 
using the KT value 77 = 1/4. The collapse of the data from 
the various system sizes is clearly excellent, and quantitatively 
confirms the KT scenario. 



scaling by 

a{L) = L-^h{LlO. XAL) - L'-"f2{L/0, (8) 

with scaling functions fi. Regarding Xo-, this implies that 
the chemical potential of the susceptibility maximum, 
see also Fig. [1] must occur at the same argument of the 
scaling function 

^ = c, (9) 



with c a constant of order unity. By substitution of Eq. ^ 
one easily derives that, to leading order, /i^ is shifted from 
|17I18|19) 



I^L — Moo 



(lnL)2 



(10) 



In Fig. [71 we have plotted versus l/(lnL)^, and the 
data are well described by Eq. IjlOp . From the fit we obtain 
A*oo = 5.187, which is also inside the region of the cumulant 
intersections of Fig. [31 In order to obtain the density poo 
at the transition, we have measured 



Pl 



(11) 



as a function of L. The result is shown in Fig.[8l where we 
have assumed that is also shifted according to Eq. (fTO|) : 
by fitting we obtain poo ~ 6.98. Note that Eq. (fTO|) proba- 
bly only approximately describes the density shift, as the 
latter is not a field variable, in contrast to the chemical 
potential. The fit, nevertheless, appears to describe the 
data well. 

We now measure the exponent rj, using the method of 
Loison [20]. The basic idea is that also the Binder cumu- 
lant is expressed by a finite size scaling form Ui — g{L/£^), 



Fig. 10. Application of Loison's method [20] to obtain the 
critical exponent /9 in a two-dimensional fluid of hard needles. 
Plotted is cr(L) Li^ versus U4,, for various system sizes L, using 
the KT value /3 = 1/8. 



with g a scaling function. Formally, this can be inverted 
^/■C = .9^^(^4); substitution into Eq.® yields 



(12) 



with h another scaling function, which could be expressed 
in terms of /2 and g, but the precise form does not matter. 
Hence, if we plot X<y{L) L^'^^^^^ versus U4, the data from 
different system sizes should collapse onto each other, pro- 
vided the correct value of is used. The result is shown in 
Fig.[9l where the KT exponent rj — 1/A was used. The col- 
lapse of the data from the various system sizes is excellent, 
and gives quantitative confirmation of the KT scenario in 
fiuids of hard needles. We observed that the quality of the 
collapse quickly deteriorates when a different exponent rj 
is used; the numerical uncertainty in r] is around ±0.01. 

Similarly, for the order parameter, we expect a data 
collapse when cr{L) versus U4 is plotted, provided the 
correct value of /? is used. The result is shown in Fig. \W\ 
where /? = 1/8 was used. Again, the collapse is very rea- 
sonable, except in the "tails" at high values of U4. Note, 
however, that here one enters the regime t < 0, where the 
scaling is expected to break down. Compared to the sus- 
ceptibility, we observed that the quality of the collapse is 
less sensitive to the precise value of (3 being used. The 
numerical uncertainty in /? is consequently larger, and 
around ±0.05. 

Finally, we discuss the compressibility Xp, see Fig. [H 
The data show the formation of a peak, which grows mildly 
with increasing system size. Note also that the maximum 
occurs well below poo of the KT transition. Interestingly, 
compared to lattice simulations of the XY model, Xp be- 
haves conform the specific heat [21|22I23|24] . For KT tran- 
sitions, the exponent of the specific heat is negative, mean- 
ing that it does not diverge. If we accept that Xp is the 
specific heat analogue, the most likely scenario is that the 
peaks in Fig. [4| saturate at finite heights in the thermody- 
namic limit. 
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Fig. 11. Plot of the order parameter "slope" dcr/d/i versus 
the chemical potential /i, for several system sizes L as indi- 
cated. Note the presence of the maximum, and also that the 
maximum increases with system size. This means that, even 
though nematic order is absent in the thermodynamic limit, 
finite-sized samples still reveal a steep rise in nematic order at 
special values of the chemical potential. 



4.2 Conventional critical scaling 

As is well known, there is no magnetization in the XY 
model in the thermodynamic limit. In finite systems, how- 
ever, the KT transition is always accompanied by a rise 
in magnetization. BH have shown that this eflFect is so 
strong, it survives in experiments |4j. Using renormaliza- 
tion group arguments, they demonstrate that the increase 
in magnetization is conform a conventional power law in 
temperature, with an associated efi^ective critical exponent 
/3eff = 37r^/128 K, 0.23. More remarkably, /?off is universal. 
Indeed, many experiments on XY-like systems yield ex- 
ponents in agreement with the BH prediction [25|- Hence, 
/3eff appears to be a genuine signature of XY universality, 
even though in the thermodynamic limit it has no mean- 
ing. 

Considering now the case of two-dimensional hard nee- 
dles, it seems reasonable to expect the validity of the BH 
scenario also. In the vicinity of the KT transition, the or- 
der parameter a rises sharply, and the slope 



A<j/Aii={{SN)~{S){N))/A, 



(13) 



attains a maximum, see Fig. [TTJ Hence, we propose to 
re-analyze our a and Xa data, but this time assuming 
conventional critical scaling. The latter is easily done in 
standard scaling plots [9]. For the susceptibility, one plots 
X(t{L) L^^'I^ versus tL^^^ , with 7 and v the critical expo- 
nents of the susceptibility and correlation lengthy respec- 
tively; for the order parameter, one plots o{V) L^'^ versus 
tL^^'^ , with P the critical exponent of the order parameter. 
Recall that t = jioo — M is the distance from the critical 
point. Provided correct values of /ioo and the critical ex- 
ponents are used, the curves from different system sizes 
are expected to collapse. 
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Fig. 12. Susceptibility scaling plot assuming a conventional 
critical point. Shown is Xo-{L) L '^^ versus tL^^" , using v ~ 
1.33, 7 « 2.33, and /ioo ~ 5.183. 
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Fig. 13. Order parameter scaling plot assuming a conven- 
tional critical point. Shown is a{L)L^^^ versus tL^^^, using 
u « 1.33, P ^ 0.18, and /^oo ~ 5.183. 



The results are shown in Fig. [12] and Fig. [131 for the 
susceptibility and order parameter, respectiveljH. The col- 
lapse looks reasonable in both cases, and we obtain v ~ 
1.33, 7 « 2.33, (3 « 0.18, and /ioo « 5.183. The estimate of 
/ioo is very close to the KT result of the previous section. 
Obviously, over the range of L available in simulations, de- 
viations from the logarithmic shift of Eq. (fTO|) over a power 
law will be small. Note also that the critical exponents ob- 
tained from the scaling plots are consistent, in the sense 
that hyperscahng j + 2(3 = dv is obeyed, with d the spatial 
dimension; substitution of our estimates yields d « 2.02. 
Hence, in agreement with BH, we find that finite systems 
of needles indeed give rise to effective critical exponents. 
However, our result /3cff ~ 0.18 seems rather far removed 
from the theoretical BH prediction /3cff,BH ~ 0.23. 



^ In these plots, the relative distance t 
actually used. 
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5 Discussion and Summary 

We have presented grand canonical simulation results of 
the IN transition of hard needles in two spatial dimen- 
sions. Our results are consistent with previous simulation 
studies of this model |1I2) , and confirm that the transition 
is of the KT type. The novelty of the present work has 
been the combination of multiple histogram reweighting 
[H] with finite size scaling. This combination facilitates an 
accurate scaling analysis, and critical exponents can be 
meaningfully obtained. Indeed, our data show that the 
XY exponents 77 = 1/4 and /? = 1/8 set in at the transi- 
tion. The chemical potential and density at the transition 
were found to be /ioo ~ 5.187 and poo ~ 6.98, which can 
be compared to Khandkar and Barma (KB) pj, who use 
deposition-evaporation dynamics to study the same tran- 
sition. In the KB case, the control parameter is the ratio 
of deposition-to-evaporation moves k; the latter is related 
to the chemical potential via k = exp(/z)//5 |2j. KB re- 
port ~ 25.8, while our estimates of ^00 and Poo yield 
Koo ~ 25.6, which is remarkably close. 

Considering the behavior of the Binder cumulant, see 
Fig.m our data reveal an approximate intersection point, 
occurring close to /ioo- This behavior is conform XY uni- 
versality |20|26l27j . Interestingly, the data of KB do not re- 
veal cumulant intersections [2]. Possibly, intersections are 
also present in the KB data but on a finer k scale, accessi- 
ble only with histogram reweighting |28^8]. Hasenbusch de- 
rived the very precise estimate lim^^oo 1/C4,kt — 1.018192 
for the value of the cumulant at a KT transition [29]; sim- 
ulation data also show that the Hmit is approached from 
above with increasing L [29J. Since our definition of the 
cumulant uses the inverse, we anticipate an increase of U4 
with increasing L. Inspection of Fig. [3] reveals a shift of 
the intersection point toward higher U4 - consistent with 
Hasenbusch - but our data still deviate from the theoret- 
ical value by w 3%. Presumably, much larger systems are 
required before the limiting value is observed, owing to 
strong subleading corrections to scaling [29] . 

We have also analyzed our data assuming a conven- 
tional critical point. The motivation was to test whether 
the effective exponent (3cS, predicted by BH, can also be 
observed. While we do recover effective exponents obeying 
hyperscahng, our estimate of /3cff does not conform to the 
BH prediction, the deviation being over 20%. Hence, the 
most consistent description of our results is provided by 
the KT scenario, in agreement with the pioneering simu- 
lations of Frenkel and Eppenga [1]. 
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